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Abstract. - We discovered an oscillatory instability in a system of inelastically colliding hard 
spheres, driven by two opposite "thermal" walls at zero gravity. The instability, predicted by a 
linear stability analysis of the equations of granular hydrodynamics, occurs when the inelasticity 
of particle collisions exceeds a critical value. Molecular dynamic simulations support the theory 
and show a stripe-shaped cluster moving back and forth in the middle of the box away from 
the driving walls. The oscillations are irregular but have a single dominating frequency that is 
close to the frequency at the instability onset, predicted from hydrodynamics. 



Introduction. - Granular gas - a system of inelastically colliding hard spheres - is a 
minimalistic model of granular flow that has received much recent attention In partic- 
ular, this model captures clustering, a generic and fascinating phenomenon in rapid granu- 
lar flows resulting from the coUisional loss of the kinetic energy of random motion of par- 
ticles. The formation and structure of granular clusters have been investigated in many 
works, both in the context of a freely "cooling" granular gas and for driven granular 
gases [3|l|3ElE||3|3[Tni[niIIlII3iniIISlIiniiniIIH|- The basic physics of clustering can often 
be described in a hydrodynamic language, in terms of a collective condensation mode driven 
by bulk losses of energy. Similar condensation processes, driven by radiative energy losses in 
gases and plasmas, have been known for a long time JHI- A complete understanding of the 
properties of granular gas (a system intrinsically far from equilibrium) is still lacking. One im- 
portant open question is the exact criteria for the validity of kinetic theory and Navier-Stokes 
hydrodynamics, developed in the 80-ies PHII^ . 

Driven granular gas has the advantage that a steady state is achievable: the energy sup- 
plied into the system can be balanced by the coUisional energy losses. We have focused in 
this work on a simple model: inelastic smooth hard spheres confined in a rectangular box and 
driven by two opposite "thermal" walls at zero gravity. Prototypical granular systems of this 
type were considered by many workers. In the early work, theoretical 5 8 10] and experimen- 
tal 0], the "stripe state": a denser and "colder" stripe-like cluster of particles, away from the 
driving wall(s), was documented. The formation of the stripe state has a simple explanation. 
Due to the inelastic particle collisions, the granular temperature goes down with an increase 
of the distance from the thermal wall(s). To maintain the steady-state pressure balance, the 
particle density should increase with this distance. For a large density contrast, the stripe 
state is observed. 
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It has been found more recently that the stripe state is only one, trivial state of this system. 
In a certain region of parameters the stripe state becomes either unstable or metastable and 
undergoes a spontaneous symmetry breaking and phase separation. As the result, a higher- 
density "phase" coexists with a lower-density phase in the direction parallel to the driving 
wall(s) |lllll2l[T^ll5[ll6[ll7lll4j . A fascinating phenomenology of this far-from-equilibrium 
phase separation, and analogy with van der Waals gas 14 put this system into a list of 
pattern-forming systems far from equilibrium. 

In this work we report a new and very different symmetry-breaking instability in the same 
system. Like the phase separation instability, the new instability occurs when the inelasticity 
q = (1 — r)/2 exceeds a critical value depending on the rest of the parameters of the system 
(here r is the coefficient of normal restitution of the particle collisions). In contrast to the phase 
separation instability, this is an instability in the direction normal to the driving wall(s), and 
it is oscillatory. As it develops, the stripe-shaped cluster in the middle of the system exhibits 
large irregular oscillations with a single dominating frequency. 

The rest of the paper is organized as follows. First, we formulate a hydrodynamic model 
and briefly describe the static stripe state. Then we present a linear hydrodynamic stability 
analysis of the stripe state that predicts oscillatory instability. The instability is then observed 
in molecular dynamic (MD) simulations. Finally, we briefly discuss our results. 

Model and static stripe state. - Let N ^ 1 identical smooth hard disks with diameter 
d and mass m move and inelastically collide inside a two-dimensional rectangular box with 
width L and height H. There is no gravity in the model. The system is driven by two thermal 
walls, with the same temperature Tq, located at a; = —L/2 and x = L/2. We assume nearly 
elastic particle collisions, g <C 1, and moderate granular densities: n/ric < 0.5, where n is 
the local number density of the particles, and ric — 2/(\/3rf^) is the hexagonal dense packing 
density. These assumptions enable us to employ a standard version of Navier-Stokes granular 
hydrodynamics {21; . Hydrodynamics is expected to be valid when the mean free path of the 
particles is much smaller than any length scale (and the mean time between two consecutive 
collisions is much smaller than any time scale) described hydrodynamically. Let us measure 

1 /2 

the distance in the units of L, the time in the units of L/Tq , the density n(r, t) in the units 
of He, the granular temperature T(r, t) in the units of Tq, and the mean velocity v(r, t) in the 
units of Tq^^ . Then the hydrodynamic equations j21ll22j become dimensionless: 



dn/dt + nV - v ^0, ndv/dt^V-P, 
ndT/dt+pV - v = eV ■ (T^/^FiVT) - eRnG T^/^ , 



(1) 



Here P = [-p + enGT^/^tr (D )] I eFaT^/^ D is the stress tensor, p = nT{l + 2G) is the 
equation of state, D — (1/2) [Vw -(- (Vu)"^] is the rate of deformation tensor, D = D — 
(1/2) tr (D ) I is the deviatoric part of D, and I is the identity tensor. The fimctions Fi and 
F2 are the following: 



. 97r / 2 
^+16^+3G 



and Ft = riG 



TT / 1 



where 
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The two scaled parameters entering Eqs. ^ are the smaU parameter e = 2TT^^^'^d/ L <C 1 and 
the parameter R ~ (16/7r)(l — r)e^^ that shows the relative role of the coUisional heat loss and 
heat conduction. The boundary conditions for the temperature become T{x = — l/2,y,t) = 
T(x = 1/2, y^t) — 1. For the velocity we demand zero normal components and slip (no 
stress) conditions at all boundaries. We assume that H is small enough, so that no symmetry 
breaking can occur in the y-direction ' 1HI12|IT^ . Therefore, the hydrodynamic variables are 
independent of y, while the mean flow, if any, is directed along the x-axis. The conservation of 

1/2 

the total number of particles yields J_i/2 '^^ n{x, t) = f , where / = N/ [LHric) is the average 
area fraction of the grains, an additional scaled parameter of the problem. 
The static stripe state is described by the equations 

Kr,(l + 2G,)]' = and (T^^ JiT,')' = i?n,G, T^/^, (2) 

where Gg — G{ns), Ji = Fi{ns), and the primes denote the a::-derivatives. Figure^ shows 
an example of the temperature and density profiles of this state. Notice that, despite the 
presence of relatively large temperature and density gradients in the stripe state, the mean 
free path of the particles remains, at q <^ I, much smaller than the characteristic length scale 
of these gradients Therefore, hydrodynamics remains valid. 




Fig. 1 - (a) - The scaled temperature profile of the static stripe state for R = 4.42-10^ and / = 0.00063. 
The inset shows the respective density profile. Except in a small vicinity of a: = 0, the gas is dilute 
in this example. As the result, the scaled temperature profile there is T ~ 2|a::|, as follows from the 
dilute limit of Eqs. lj2j. (b) - The solid line is the oscillatory instability threshold R = Rc versus 
the average area fraction / at e = 4 ■ 10~^. The instability occurs when R > Rc- The dashed lines 
show the borders of the spinodal interval of the phase separation instability that develops in the same 
system at large enough H. The two asterisks correspond to two MD simulations (see Figs. |2|-0J. 

Linear stability analysis. ~ The linear stability analysis involves linearization of Eqs. (|Hl 
around the static profiles ns{x) and Ts{x). We substitute 

n{x,t) ^n,{x) +e-'''^^iy{x), T{x,t) = T^ix) + e~''^* Q{x), and v{x,t) = e'"'^* ^{x), (3) 

where Imw > (< 0) corresponds to growth (damping) of the small perturbation, and obtain 
the linearized equations: 



— iujv + (risipy — , 
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lojn,^ = ~ (n,(l + 2 Gs) Q)' ~ (T, J3 i^Y + e (t^^ V^')' + | {t^'' J2 V'')' , 
icj e + Hs %l> T'^ + T, (1 + 2 G,.) ^' 

(Ty2 e)' jj ' + e f Tl'^ T',J,A' -end Tl'^ G, e + J5 , (4) 



V 

where J2 = -F2("-s), >/3 = 1 + 2Gs + 2nsdGs/dns, 

■h = [4^ + 127rG, + (16 + 9^)G2] + ^ [-47r + (16 + 97r)G2] ^ , 

>/5 = Gs + nsdOs/dus and dGs/dug — (97r/4)(16-\/3 + 7rns)(6 — VStttt-s ) • Eliminating i^(x) 
from equations Q , we arrive at a linear eigenvalue problem for the two-component eigenvector 
U(a;) — [^'(2:), 0(2;)], corresponding to the complex eigenvalue uj: 

AU" + BU' + CU = 0. (5) 

The elements of matrices A, B and C are 

All = -{ns/iuj)Ts J3 + eaofisGs + {e/2) ao J2, A21 = e {ns/iu!)T^^^ T!, J4, 

A12 = 0, A22 = eao Ji, 

Bii = -{ns/iuj) Tl J3 - (us/iuj) Tg n'^ Jq - {2n'Jiuj) Tg J3 + eoi G^ + eag n'^ J5 

+ (e/2)ai J2 + (e/2)ao< J7, 

B21 = ~nsTs{l + 2Gs)+e{ns/iu){Tl'^T'g)' .h + e{ng/iuj)Tl'^T'gn'gJs 

+ e {2n'Jiuo) T}'^ J4 - e R {us/iuj) T^/^ J5, 

B12 = -ns{l + 2Gs), B22 = 2faiJi + eaon'g Ji, 

Gil = iujris- {n'Jiuj)T^J3- {{n'^)^/iuj)TsJ6~ {n'^/iuj)TsJ3, 

C21 = -nsT'g + e{n'Jiu;){T}/^T'JJ, + ein'JiLo)T}/^T'gn'gJs 

+ € (n'l/iuj) Tl'^ r; Ji-eR [n'Jiuo) T^'^ J5, 

G12 = -n'gJz, G22 = «wns + £02 Ji + eai n(, J4 - (3/2) ei?ns Gs T^^, (6) 



and 



Je = AdGs/dus + 2nsd'^Gs/dTil, 

Jj = {l/8)G-^[TT + 2TTGs + {8 + Tr)G'^,]+{l/8)G-^ns[-TT + {8 + Tr)G'^g]dGs/dns, 

Js = {^/2)G-^ng{dGs/dnsf -{TT/^)G-^{2dGs/d7i,+ngd^Gs/d7il) 

+ (1/16)(16 + 97r) (2 dGs/dus + d^Gs/dnl), 

d^Gg 



dnl 



= (97rV2) (75 + %/3 7rn,) (6 - VS 
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We have denoted for brevity oq = Ts (x) , ai ~ u'q and 02 = «[,'. The boundary conditions 
for the functions O and ip are the following: 

e(±l/2) = V(±l/2)-0. (7) 

Equations ©-0 define a linear boundary-value problem: there are two boundary conditions 
at one wall, and two at the other wall. A simple and accurate numerical algorithm, that we 
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Fig. 2 - A snapshot of the system of A'' = 3390 inelastic hard disks within the instability region, 
obtained in a MD simulation. The left and right walls are thermal walls. The coefficient of normal 
restitution r = 0.85. See the text for the rest of the parameters. A large deviation of the cluster from 
the center of the system a; = is clearly seen. 



developed earlier enables one to avoid the unpleasant shooting in two parameters. The 
algorithm yields the complex value of lu. Varying R at fixed / and e, we determined the 
critical value R — Rc for the instability onset from the condition Imcj = 0. The instability 
occurs at R > Rc- Given the rest of the parameters, the instability occurs when q exceed a 
critical value. Importantly, the instability is oscillatory, as Heuj at the onset. Figure ^3 
shows the critical value R = Rc versus the average area fraction / at a fixed e. 

MD simulations. ~ To verify the hydrodynamic predictions, and to follow the dynamics 
of an oscillating stripe in a nonlinear regime, we performed MD simulations with N = 3390 
inelastic hard disks with unit mass and diameter. As Fig. ^ shows, the oscillatory instability 
region lies within the spinodal interval of the phase separation instability |llll3ll4j . The phase 
separation instability is suppressed when the aspect ratio HjL is less than some threshold 
value depending on R ^1L..12..13) . It has been found recently that the system exhibits large 
fluctuations even well below the threshold jJHI- Therefore, to isolate the oscillatory instability 
from the phase separation instability, we worked with a very small aspect ratio: the box 
dimensions were R — 208 and L — 2% 209. The parameters / = 5 • 10~^ and e = 4 • 10~^ were 
fixed. The parameter R varied from 4.42 • 10^ to 9.55 • 10^. This was achieved by varying r from 
0.85 to 0.997. We employed the event-driven algorithm described in the book of Rapaport 
Upon collision with a thermal wall, the normal component of the particle velocity was drawn 
from a Maxwell distribution with Tq = 1, while the tangential component of the velocity 
remained unchanged. The initial spatial distribution of the particles was uniform; the initial 
velocity distribution was Maxwell's with Tq = 1. Figure El shows a typical snapshot of the 




0.2 0.4 0.6 J 0.8 1 0.2 0.4 0.6 ^ 0.8 1.0 

Fig. 3 - The time dependence of the center-of-mass coordinate Xcm , obtained in MD simulations with 
A'^ = 3390 particles within (a) stable region, r — 0.997 and (b) unstable region, r = 0.85. The time 
unit is 10^ tMD, where tMD = d/Tg^^ — 1. 
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Fig. 4 - The power spectra within the stable (a) and unstable (b) regions. The frequency unit is 
2-K ■ 10~^. The parameters in Figs, a and b are the same as in Figs. andj^a, respectively. Notice 
the large differences in scales. 

system obtained in a MD simulation within the instability region. A large deviation of the 
cluster from the center of the system a; = is clearly seen. A movie of this simulation can be 
seen at http:/ /huji-phys. phys.huji.ac.il/staff/Khain/abstract. html The oscillatory motion of 
the cluster was monitored by following the x-coniponent of the center of mass of all particles, 
Xcm, versus time. Figure 13 shows this dependence as observed in the stable (a) and unstable 
(b) regions. The large-amplitude low-frequency irregular oscillations in the unstable case 
are clearly distinguishable from the small-amplitude high-frequency noise observed in the 
stable region. The power spectra for these two cases are shown in Fig. 01 Notice the large 
difference in scale in the horizontal axes and the huge difference in the vertical axes. In 
contrast to the low-amplitude noise of Fig. the power spectrum within the instability 
region (Fig. ^p) has a sharp peak at a single frequency 1.03 • 10""*. This frequency is fairly 
close (within 30%) to the hydrodynamic frequency at the instability onset which, in the units 
of f^j-) = 1, is 1.45 • 10~^. These two frequencies are not expected to coincide. Firstly, the 
steady-state oscillations, observed in the MD simulation (Figs, and^), must be affected 
by hydrodynamic nonlinearities that are neglected in the linear stability analysis. Secondly, 
the next-order corrections in q can become significant at the moderate value of r = 0.85 used 
in this simulation |22) . 

The presence of a significant continuum part in the power spectrum of Xcm{i) in Fig- Eb 
clearly shows that the cluster oscillations are chaotic. The same conclusion follows from the 
analysis of the phase trajectory in the phase plane {Xcm, -^cm)- The chaotic component of the 
oscillations is due to irregular cluster motion, not due to the dynamics of the surrounding gas. 
We double-checked it by following the position X^, of the maximum number density of the 
system (integrated over the y-direction) versus time. The graph oiX^(t) almost coincides with 
the graph of Xcm{t)', the respective power spectra are almost indistinguishable. Obviously, 
the chaotic component of the cluster oscillations is completely missed by the linear stability 
analysis. Furthermore, hydrodynamics predicts, at fixed / and e, a sharp instability onset 
at R ^ Re- On the contrary, a fairly smooth transition is observed in MD simulations. We 
believe it is due to the relatively small number of particles that causes significant fluctuations 
and smoothes the transition [221 ■ These fluctuations make it difficult to verify the reentrant 
behavior (see the right branch of the solid line in Fig. ^p) predicted by hydrodynamics. 

Summary and Discussion. - We have discovered a new oscillatory instability in a system 
of inelastically colliding hard spheres, driven by two opposite "thermal" walls at zero gravity. 
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The instability has a hydrodynamic character and can therefore be interpreted in terms of the 
collective modes of the system. There are three such modes, coupled due to the inhomogeneity 
of the unperturbed state of the system. Two of the modes are oscillatory acoustic modes, the 
third one is a non-oscillatory entropy mode (the fourth, shear mode, is suppressed because 
of the very small aspect ratio of the box). It is the acoustic modes that become unstable at 
R > Rc- To some extent, this instability is similar to the acoustic instability (the so called 
"overstability") found in externally heated and radiatively cooling optically thin plasmas 'TS]. 
The future work should address the nonlinear regime of the cluster oscillations and clarify 
the role of discrete-particle noise versus the spatio-temporal chaos caused by hydrodynamic 
nonlinearities. 

* * * 
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